```{r}

rm(list = ls())

library(pacman)
p_load("fixest", "haven", "conleyreg", "readxl", "dplyr", "lmtest", "ggplot2", "ggthemes", "jtools")

remotes::install_github("shuo-zhang-ucsb/did_multiplegt")
Panel_long <- read_xlsx("Panel_28Jan2024.xlsx")
Panel_long <- read_xlsx("\\\\tsclient/C/Users/tanay/Dropbox/ADS/Panel_28Jan2024.xlsx")

```


##DID_MGT Regressions
```{r}

Panel_long$First_treatment <- ifelse(is.na(Panel_long$First_treatment), 9999,Panel_long$First_treatment)

Panel_long$ADS2 <- ifelse(Panel_long$year == Panel_long$First_treatment - 2 | Panel_long$year == Panel_long$First_treatment - 1, 1, Panel_long$ADS)

HDI_1 <- did_multiplegt(Panel_long, Y = "HDI_scaled", G = "code", "year", D = "ADS2", controls = c("Pop", "dist", "neigh"), cluster = "code", brep = 1, dynamic = 16, covariance = T)

EDI_1 <- did_multiplegt(Panel_long, Y = "EDI_scaled", G = "code", T = "year", D = "ADS2", controls = c("Pop", "dist", "neigh"), cluster = "code", brep = 500, dynamic = 16, covariance = T)

HDI_1_p <- did_multiplegt(Panel_long, Y = "HDI_scaled", G = "code", T = "year", D = "ADS", controls = c("Pop", "dist", "neigh"), cluster = "code", brep = 500, placebo = 14, parallel = T)

EDI_1_p <- did_multiplegt(Panel_long, Y = "EDI_scaled", G = "code", T = "year", D = "ADS", controls = c("Pop", "dist", "neigh"), cluster = "code", brep = 500, placebo = 14, parallel = T)


```


##EDI by type
```{r}

EDI_accident <- did_multiplegt(Panel_long, Y = "EDI_accident_scaled", G = "code", T = "year", D = "ADS2", controls = c("Pop", "dist", "neigh"), cluster = "code", brep = 500, dynamic = 16, covariance = T,parallel = T)

EDI_accident_p <- did_multiplegt(Panel_long, Y = "EDI_accident_scaled", G = "code", T = "year", D = "ADS", controls = c("Pop", "dist", "neigh"), cluster = "code", brep = 500, placebo = 14,parallel = T)

EDI_other <- did_multiplegt(Panel_long, Y = "EDI_other_scaled", G = "code", T = "year", D = "ADS2", controls = c("Pop", "dist", "neigh"), cluster = "code", brep = 500, dynamic = 16, covariance = T,parallel = T)

EDI_other_p <- did_multiplegt(Panel_long, Y = "EDI_other_scaled", G = "code", T = "year", D = "ADS", controls = c("Pop", "dist", "neigh"), cluster = "code", brep = 500, placebo = 14,parallel = T)


```


##Plotting the DID_MGT results
```{r}
p_load("broom")

tidy <- function(x, level = 0.95) {
  ests = x[grepl("^placebo_|^effect|^dynamic_", names(x))]
  ret = data.frame(
    term      = names(ests),
    estimate  = as.numeric(ests),
    std.error = as.numeric(x[grepl("^se_placebo|^se_effect|^se_dynamic", names(x))]),
    N         = as.numeric(x[grepl("^N_placebo|^N_effect|^N_dynamic", names(x))])
    ) |>
    # For CIs we assume standard normal distribution
    within({
      conf.low  = estimate - std.error*(qnorm(1-(1-level)/2))
      conf.high = estimate + std.error*(qnorm(1-(1-level)/2))
      })
  return(ret)
}


Regs <- list(
HD = HDI_1, HD_p = HDI_1_p, ED = EDI_1, ED_p = EDI_1_p, ED_accident = EDI_accident, ED_accident_p = EDI_accident_p, ED_other = EDI_other, ED_other_p = EDI_other_p)


Regs <- sapply(Regs, tidy, simplify = F, USE.NAMES = T)

# for (i in 1:length(Regs)){
#  filename = paste0(names(Regs)[i], ".xlsx")
#  writexl::write_xlsx(regs[[i]], filename)
# }

Plots <- list()
for (i in 1:length(Regs)){
  
  
  xlabl <- ifelse(substr(names(Regs)[i], start = nchar(names(Regs)[i]), stop = nchar(names(Regs)[i])) == "p", "Years since implementation of ADS -2", "Years since implementation of ADS -3") 
  
  ylabl <- ifelse(substr(names(Regs)[i], start = 1, stop  = 2) == "HD", "Effect of ADS on human deaths", "Effect of ADS on elephant deaths")
  
  
  Plots[[i]] <-  Regs[[i]] |>
  within({
    term = gsub("^placebo_", "-", term)
    term = gsub("^effect", "0", term)
    term = gsub("^dynamic_", "", term)
    term = as.integer(term)
    }) |>
    ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() + geom_hline(yintercept = 0) + labs(x = xlabl, y = ylabl) + theme_clean(base_size = 16) + theme(axis.text.x=element_text(size=18), axis.text.y=element_text(size=18), axis.title.x = element_text(size=18), axis.title.y = element_text(size=18))


  ggsave(filename = paste0("./", names(Regs)[i], ".jpg"), plot = Plots[[i]], device = "jpeg")
}

```


##function for calculating average effects
```{r}

Avg_effect <- function(x, level = 0.95) {
  ##retrieve the covariances (17C2 = 17*8 terms)
  ests = x[grepl("^cov_effect", names(x))]
  covars = data.frame(
    term      = names(ests),
    covar     = as.numeric(x[grepl("^cov_effect", names(x))])
    )
    ##drop first three years (these are the first 17C2 - 14C2 = 45 terms)
    covars <- covars[46:nrow(covars), ]
  ##retrieve the standard errors and calculate the variances
  ests2 = x[grepl("^se_", names(x))]
  vars = data.frame(
    term      = names(ests2),
    se        = as.numeric(x[grepl("^se_", names(x))])
    ) |> within({var = (se^2)})
    ##drop first three years
    vars <- vars[4:nrow(vars), ]
  ##retrieve the coefficients
  ests3 = x[grepl("^effect|^dynamic_" , names(x))]
  effects = data.frame(
    term      = names(ests3),
    effects   = as.numeric(x[grepl("^effect|^dynamic_", names(x))])
    )
    ##drop first three years
    effects <- effects[4:nrow(effects), ]
  ##retrieve the number of observations
  ests4 = x[grepl("^N_dynamic|^N_effect" , names(x))]
  nobs = data.frame(
    term      = names(ests4),
    nobs   = as.numeric(x[grepl("^N_dynamic|^N_effect", names(x))])
    )
    ##drop first three years
    nobs <- nobs[4:nrow(nobs), ]
  ##calculate the cross products of observations
    N1N2 <- c()
    for (i in 1:nrow(nobs)){
      if (nrow(nobs)-i > 0) {
        for(j in 1:(nrow(nobs)-i))
          N1N2 <- append(N1N2, nobs[i, 2]*nobs[i+j, 2])
  }
    }
  ##calculate the average effect for t = 3 to t = 17 and its se
    effects$nobs <- nobs$nobs
    effects <- mutate(effects, abeta = effects*nobs)

    covars$ab <- N1N2
    covars <- mutate(covars, abcovar = ab*covar)

    vars$nobs <- nobs$nobs
    vars <- mutate(vars, N2var = (nobs^2)*var)

    avg_effect = data.frame(
    coef = sum(effects$abeta)/sum(nobs$nobs),
    se = sqrt(sum(vars$N2var) + 2*sum(covars$abcovar))/sum(nobs$nobs))
    avg_effect <- mutate(avg_effect,
                           ci_low =  coef - se*(qnorm(1-(1-level)/2)),
                           ci_high = coef + se*(qnorm(1-(1-level)/2)))
  return(avg_effect)
}

Avg_effect_HD <- Avg_effect(HDI_1)
Avg_effect_ED <- Avg_effect(EDI_1)

Average_effects <- tibble(model = c("Human deaths", "Elephant deaths"), rbind(Avg_effect_HD, Avg_effect_ED))

DIDMGT_plot <- ggplot(data = Average_effects, aes(x = model, y = coef)) + geom_point() + geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.3) + geom_hline(aes(yintercept=0), linetype="dashed") + xlab("Model") + ylab("Average effect of ADS") + theme_clean()

ggsave(DIDMGT_plot, filename = "./DID_MGT_average_effects.jpeg", scale = 4)


```

